Discreteness-Induced Criticality in Random Catalytic Reaction Networks 
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Universal intermittent dynamics in a random catalytic reaction network, induced by smallness in 
the molecule number is reported. Stochastic simulations for a random catalytic reaction network 
subject to a flow of chemicals show that the system undergoes a transition from a stationary to 
an intermittent reaction phase when the flow rate is decreased. In the intermittent reaction phase, 
two temporal regimes with active and halted reactions alternate. The number frequency of reaction 
events at each active regime and its duration time are shown to obey a universal power laws with the 
exponents 4/3 and 3/2, respectively. These power laws are explained by a one-dimensional random 
walk representation of the number of catalytically active chemicals. Possible relevance of the result 
to intra-cellular reaction dynamics is also discussed. 
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Most intra-cellular reactions progress with the aid of 
catalysts. These reactions are essential for a cell to sur- 
vive and grow. All catalysts that are proteins are syn- 
thesized within a cell as a result of catalytic reactions. 
Thus, researchers have started investigating catalytic re- 
action networks in order to develop a theoretical model 
of intra-cellular dynamics andprotocells for their role in 
the origin of life[!|, i, % i, i, i, Si, 0; research in this 
field was pioneered by KauffmanQ. Recent studies on a 
growing cell model comprising random catalytic reaction 
networks have also revealed universal statistical laws on 
chemical abundances. These results are confirmed by the 
gene expression data of the present cells[^|3|- 

Generally, cells consist of a large number of chemical 
species, some of which flow in and out through the mem- 
brane. There are also other chemical species that are 
not present in a large quantity. However, some chem- 
ical species play an important role even at extremely 
low concentrations amounting to only a few molecules 
per cell (IT 
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In such cases, the fluctuation 
and discreteness in the molecule number are very im- 
portant as far as temporal variations and steady distri- 
bution of chemical concentrations are considered. In- 
deed, recent studies on catalytic reaction dynamics have 
demonstrated that discreteness in molecule number can 
induce drastic changes in the temporal variations and the 
steady distribution of chemical concentrations obtained 
from the rate equation model 14, 1^, 1^, 17 1. However, 
in most thermodynamics studies, the molecule number is 
assumed to be large while the number of species is as- 
sumed to be rather small. In order to study the catalytic 
reaction network of a cell, it is important to consider 
the case a large number of chemical species and a small 
molecule number. 

In this study, we investigate a simple model of a ran- 
dom catalytic reaction network subjected to a flow of 
chemicals. In a region with a weak flow, where the dis- 



(a) 



(b) 



Inflow ■ 



9. 



^ f <^ 3/ X n cjj 



•® 10 



FIG. 1: (a) An illustration of a catalytic reaction B -\- C - 
A-\- C, and (b) an example of a catalytic reaction network. 



creteness in molecule number is large, it was observed 
that the reaction events occurred intermittently, sepa- 
rated by a quiescent state. It was also observed that 
the number distribution of reaction events during each 
reaction state obeyed a universal power law, display- 
ing a feature of criticality induced by the discreteness 
in the molecule number. By referring to the study of 
self-organized criticality, as investigated in models for 
sand-piles, earth qua kes, interface depinning in random 
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23|, we will explain 
the origin of this universal power law behavior. Finally, 
we discuss in brief the possible relevance of our result to 
the power laws of the distribution of residence time in 
the quasi-stationary state observed in behaviors of cells. 

We introduce a network of elementary, two-body 
catalytic reactions among a large number of chemical 
species [Q]. Assuming that chemicals in the system are 
well stirred, we discard the spatial dependence of concen- 
tration so that the state of the system can be represented 
by a set of numbers (ni, 712, ....riM), where ni(0, 1, ...) in- 
dicates the number of molecules of each chemical species i 
(1 < J < M), with M being the total number of chemical 
species. Each chemical is transformed to another by (sev- 
eral) reactions. These reactions are catalyzed by another 
chemical, i.e., the reaction from chemical B to chemical 
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A is catalyzed by a third chemical C (see Fig. 1). For 
simplicity, the reaction rates r of each reaction are set to 
be identical (the inhomogeneous reaction rate case is dis- 
carded as it does not affect the result to be discussed). 
The growth rate of molecule number ua (or the decay 
rate of molecule number ns) is given by rnsnc /V , on 
the average, where V represents the volume of the sys- 
tem. 

The entire catalytic reaction network consists of re- 
actions whose reaction paths are chosen randomly (and 
then fixed). The average number of reaction paths from 
a chemical i catalyzed by a chemical j is set to a given 
connection number K. In this study, we only consider the 
case where Kc < K << M; Kc indicates the critical con- 
nection number at the percolation threshold in random 
networks {Kc ~ log(M)) We have not considered an 
auto-catalytic reaction of the form B + C ^ 2C because 
usually such a reaction is not elementary and is realized 
as a result of a series of (non-auto-catalytic) elementary 
reactions. 

Next, we consider the flow of chemicals into a system 
from a molecular bath. We assume that Mi„ chemical 
species can flow into the system. At each instance, a 
single molecule from one of the Mi„ chemical species is 
added to the system in the rate Q/Min, where Q indi- 
cates the inflow rate of molecules. For simplicity, the 
flow rates of all chemical species are assumed to be the 
same; however, this assumption can be relaxed. We also 
assume that Mout chemical species flow out (or are de- 
composed) from the system. Because of this inflow and 
outflow of chemical species, the total number of molecules 
N = Uj varies with time, which is in contrast to our 
previous study p/7,]. 

In this study, we assume that Mi„ is sufficiently larger 
than ^JM/K. It should be noted that the average num- 
ber of reaction paths through which the Mi„ 'input' 
chemicals are catalyzed by one of themselves is given as 
MinK X Min/M. This number is much greater than unity 
when Min » a/ M/K. As long as this condition is sat- 
isfied (empirically, Min ^ '^\^M/K), a chemical reaction 
takes place even if the total molecule number is null ini- 
tially. With this constraint on Min, the behavior of the 
system is observed independent of Min ■ Here, we mainly 
show the results of the system behavior for A/i„ — M. 
Our second assumption is that Mout is much smaller than 
M. If the value of Mout is comparable to that of M, the 
number of extinct chemical species increases. As a result, 
the number of surviving chemical species is much smaller 
than M, thereby reducing the effective system size con- 
siderably. On the other hand, if this assumption is true, 
the qualitative behavior of the system is determined by 
M. In this study, we mainly present the results of the 
system behavior for Mout — 1- 

We perform stochastic particle simulation in order to 
study the possible effects of stochasticity and discrete- 




FIG. 2: (Color online) Typical temporal evolutions of N{t) 
(blue) and RN{t) (red) for M = 300 and K ^ 12 with (a) 
g = 0.3 and (b) Q = 0.001. 

ness in molecule numbers. The simulation procedure is 
as follows: At each simulation step, we randomly select 
a pair of molecules. If the pair consists of a substrate 
and one of its catalysts, then, according to the rule of 
the catalytic reaction network, the substrate molecule 
is replaced by the product molecule with the proba- 
bility r. The molecule number rii of input chemical 
(1 < i < Afin) is incremented by one with a probability of 
Q/Min per unit time. When the output chemical species 
j {{M — Mout) < j ^ M) is generated by the reaction or 
the inflow of molecules, it is removed immediately. We 
have adopted Gillespie's direct method in order to carry 
out numerical simulations for this process[25]. 

Here, we present the results of numerical simulations 
obtained from a variety of random networks and several 
sets of parameters. Without losing generality (by redefin- 
ing the time unit), we can assume r = 1 and V — 1. Fig- 
ure 2 shows the plots of typical temporal evolutions of the 
total number of molecules in the system, N{t), and the 
reaction number, RN{t), for i = 0, 1, 2, ... for M = 300, 
and K =12 with (a) Q = 0.3 and (b) Q = 0.001. Here, 
RN{t) is defined as the number of reaction events that 
have occurred between t and < + 1. When Q is large, 
N and RN show stationary, relatively small fluctuation 
around their mean values. On the other hand, when Q 
is small, there are intermittent bursts in RN , while N 
increases gradually, until the increase is replaced by the 
drastic drop in iV as a result of such bursting reactions. 
These two regimes alternate. 

The reason for change in behaviors of N{t) and RN{t) 
because of the decrease in Q is as follows. When N is 
not large enough due to low inflow rate Q, the molecules 
therein do not react with each other because catalyst 
molecules for each of existing chemical species are absent. 
The reaction then stops and RN is constant at 0, until a 
flow of necessary catalysts restarts the reaction process. 
When the reactions start, the process of production and 
elimination of catalysts is repeated until the catalysts for 
all molecules in the system disappear. 

If Q is sufficiently large, the system is in a steady state 
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FIG. 3: (Color online) p, k, and function of Q for 

several M, K, Min, and Mout- p, and k are computed by 
sampling with a step size of 10^. 



expected from the fixed point solution of the continuous 
rate equation. There exists a deviation in N and RN by 
stochastic inflows and the reaction process to produce the 
chemicals that escape from the system. However, these 
deviations are not substantial. Thus, a steady state with 
large TV is sustained so that the reaction does not stop. 

We have plotted N, the long time average of N ver- 
sus the change in Q for statistical characterization of 
the transition mentioned above. Figs. 3 (a) and (b) 
are the log-log plots of p = N'^KMout/M'^ versus Q, 
whereas Figs. 3 (c) and (d) are the plots of k = and 
k' = K — 1 as functions of Q for several M, K, = M 
and ^ M/K, and Mout = 1, 3, and 6. Here, Figs. 3(a) 
and (c) show the results for Min — M and Mout — 1 and 
Figs. 3(b) and (d) show the results under M = 1000. 

In Figs. 3(a) and (b), p « Q holds for large Q ^ 
Mout/K. This can be explained as follows. In the steady 
state, the average molecule number for each chemical is 
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FIG. 4: (Color online) Frequencies of S for M = 1000, K = 19, 
and 36, Q = 0.0003 with (a) Mi„ = M and Mout = 1 and 
(b) Min ~ M/K and Mout = 3 (M,„ = 55 for A' = 18 and 
Min ~ 27 for K — 36). These histograms are computed by 
sampling with a step size of 10*. 



N/M, which also gives the number of catalysts per reac- 
tion path. Then, the rate of production of each chemical 
is ^ K{N/M)'^. Thus the outflow rate is estimated as 
Mout X K{N /MY = p which should balance with the 
inflow rate Q at the steady state. 

On the other hand, for a small Q, N is constant 
independent of Q. This is explained by the estimate 
of the number of molecules for the discreteness-induced 
transition 17|. Here, the reactions tend to stop if the av- 
erage number of reaction paths that have corresponding 
non-zero catalyst for each chemical is smaller than unity, 
i.e., if 1. In this case, the outflow of molecules 



is suppressed; however, N increases because of the in- 
flow of molecules, until it satisfies the condition ~ 1. 
Hence, N remains at a constant value that is slightly 
larger than M/K for small Q. Indeed, Figs. 3(c) and (d) 
for small Q indicate that the Q versus k and Q versus 
k' relations for a large K are fitted by a single curve for 
each M and Mout, independent of K and Min, while k! 
approaches with an increase in M and Mout- More- 
over, the crossover point between the above two states is 
estimated as Qc Mout/K. 

We now focus on the statistical properties of the sys- 
tem at the regime of small Q with a constant N . In 
this regime, the inflow of molecules continues even af- 
ter the successive reactions have been completed. This 
situation will continue until an inflow triggers the next 
successive reactions, as shown in Fig. 2(b). Then, we 
focus on statistical properties of the reaction size S- S 
is defined as the sum of RN{t) between the interval of 
quiescent states where the molecules therein can not re- 
act with each other. Figure 4 shows the frequency of the 
reaction size S, for a system with M — 1000, K — 18 
and 36, Q = 0.0003, with M,„ = M and Mout = 1 for 
(a) and Min ~ M/K and Mout — 3 for (b). As shown, 
the frequency distribution of S, P{S), obeys a univer- 
sal power law P{S) ^ with the exponent 7^4/3 
for small enough Q- The same universal power law is 
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FIG. 5; (Color online) (a) Frequencies of D for M = 1000, 
Q = 0.0003, A'/i„ = M, and Mout = 1 with K = 1^ and 
/iT = 36 and (b) D — S characteristics for the same conditions. 
The histogram (a) is computed by the sampling with a step 
size of over 10* and (b) is computed with a step size of 2 x 10^ . 

satisfied even for other values of M, Mi„, Mo„t, and K. 

Now we discuss the origin of this universal power law 
for a small Q. Let us split the chemical species into active 
and inactive ones. An active chemical species is defined 
as that species which is present in the system along with 
its catalysts. When MA{t) = 0, the system is said to be 
in the inactive state, and when Myi(t) > 0, the system is 
said to be in the active state. 

During the inactive state, MA{i) increase due to the in- 
flow of some molecules. Reactions start taking place due 
to this inflow and production and elimination of several 
chemicals (catalysts) are repeated leading to an increase 
and then a decrease in M^(t). This process continues 
until MA{t) returns to 0, i.e., when the catalysts for all 
the molecules in the system have disappeared. Let us de- 
fine the duration of reactive state D as the time interval 
with MA{t) > 0. Note that RN{t) is considered to be 
proportional to MA{t) on average. 

Since the system under consideration is a random cat- 
alytic reaction network, we assume that the change of 
Ma is approximated by a one-dimensional random walk 
with < Ma < M. Then, the frequency of D is given by 
P{D) - Z?-" = on the basis of the first return 

time distribution for the one-dimensional random walk. 
The reaction size S, which is proportional to the area 
below the random walk before its first return, scales as 
S DP ^ 1)3/2. rpj^g^^ ^-^^ exponent 7 for P{S) - S"'^ , 
is given by the standard relation between the exponents 
a, f3, and 7: 7 = 1 -h (a - l)//3 = 4/3[2l|. 

In order to confirm the validity of the above argument, 
we have examined the frequency of D and the relation- 
ship between D and S. Figure 5 (a) shows the frequency 
of D and Fig. 5 (b) shows the D — S characteristics 
for K = 18 and K ^ 36 with M = 1000, M,„ = M, 
Mout = 1, and Q = 0.0003, which are obtained from the 
same simulation as the data plotted in Fig. 4(a). The 
Figs. 5 (a) and (b) support the relation P{D) ~ D"^/^ 
and S ~ D^/^. Note that a similar random- walk descrip- 



tion was adopted in a model of self-organized criticality 
in the anisotropic interface depinning in a quenched ran- 
dom medium [2l|. 

In this study, the dynamic aspects of a random cat- 
alytic reaction network subjected to a flow of chemi- 
cals was studied. With a decrease in the inflow rate, 
the system undergoes a transition from a stationary to 
an intermittent reaction state; however, discreteness in 
molecule number in the latter state is essential. The 
frequency of the reaction size and duration of the re- 
active states obey the universal power laws with expo- 
nents 4/3 and 3/2, respectively. Note that this critical 
behavior is obtained without any tuning parameter val- 
ues, and as long as the inflow rate is small enough to 
maintain the discreteness (0,1..) in the molecule num- 
ber. In other words, discreteness induced self-organized 
criticality (SOC) [11, [3, H HE [ll] . 

Power law of the residence time distributions in a cer- 
tain state are often reported in the behavior of cells, 
such as the residence time at a counter-clockwise ro- 
tations of flagella motor in E-coli[2^, the frequency of 
the inter-event intervals for spontaneous vesicular release 
at a Xenopus neuromuscular junction [2^. and the fre- 
quency of the inter-beat time interval in cardiac muscle 
cellj28:], where non-stationary behaviors are a result of 
intra-cellular reaction dynamics. Since our random re- 
action network model is rather simple, we cannot make 
a detailed comparison with these experiments; however, 
the discovery of a universal power law may help us to 
understand the ubiquity of non-stationary intra-cellular 
dynamics and power law behaviors. 

This research of A. A. was supported in part by 
a Grant-in Aid for Young Scientist (B) (Grant No. 
19740260). 
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